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Abstract. In the first quarter of 2006 Chicago Board Options Exchange (CBOE) introduced, as 
one of the listed products, options on its implied volatility index (VIX). This created the challenge 
of developing a pricing framework that can simultaneously handle European options, forward-starts, 
options on the realized variance and options on the VIX. In this paper we propose a new approach 
to this problem using spectral methods. We use a regime switching model with jumps and local 
0^ ' volatility defined in [T] and calibrate it to the European options on the S&P 500 for a broad range 

. of strikes and maturities. The main idea of this paper is to "lift" (i.e. extend) the generator of the 

, underlying process to keep track of the relevant path information, namely the realized variance. The 

^-^1 lifted generator is too large a matrix to be diagonalized numerically. We overcome this difficulty 

^2 ' by applying a new semi-analytic algorithm for block-diagonalization. This method enables us to 

evaluate numerically the joint distribution between the underlying stock price and the realized 
CO . variance, which in turn gives us a way of pricing consistently European options, general accrued 

variance payoffs and forward-starting and VIX options. 

d 

^ CT*. In recent years there has been much interest in trading derivative products whose underlying is 

a reaUzed variance of some hquid financial instrument (e.g. S&P 500) over the life of the contract. 
The most popular payoff functionj^ are linear, leading to variance swaps, square root, yielding 
volatility swaps, and the usual put and call payoffs defining variance swaptions. 
' It is clear that the plethora of possible derivatives on the realized variance is closely related 

1^ I to the standard volatility-sensitive instruments like vanilla options, which are also exposed to 
I other market risks, and the forward-starting options which are almost pure vega bets and are 
, mainly exposed to the movements of the forward smile. Recently Chicago Board Options Exchange 



1. Introduction 



(CBOE) introduced options on the volatility inde^ij (VIX) which are also important predictors 
^ . for the future behaviour of implied volatility. The main purpose of this paper is to introduce 
. a framework in which all of the above financial instruments (i.e. the derivatives on the realized 
variance as well as the instruments depending on the implied volatility) can be priced and hedged 
consistently and efficiently. 

Our central idea is very simple and can be described as follows. We use a model for the underlying 
that includes local volatility, stochastic volatility (i.e. regime switching) and jumps and can be 
calibrated to the implied volatility surface for a wide variety of strikes and maturities (for the case 
of European options on the S&P 500 see Figured]). The stochastic process for the underlying 
is a continuous-time Markov chain, as defined in [1], where it was used for modelling the foreign 
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^For the precise definition of these products see Appendices IB. II and IB. 2l 

^For a brief description of the these securities see Appendix IB. 31 For the definition of VIX see [12] . 
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exchange rates. The underlying process is stationary as can be inferred from the fact that the 
impHed forward volatihty smile behaves in a consistent way (see Figures [71 [8] [9l and I lOp . 

There are two features of this model that make it possible to obtain the distributions of the 
future behaviour of the implied volatility and the realized variance of the underlying. The first 
feature is the complete numerical solubility of the model. In other words spectral theory provides a 
simple and efficient algorithm (see [1], Section 3) for obtaining a conditional probability distribution 
function for the underlying between any given pair of times in the future. This property is sufficient 
to determine completely the forward volatility smile and the distribution of VIX for any maturity. 
The second feature of the framework is that it allows us to describe the realized variance of the 
underlying risky security as a continuous-time Markov chain. This makes it possible to deal with 
the realized variance by extending the generator of the original process and thus obtaining a new 
continuous-time Markov chain which keeps track simultaneously of the realized variance and the 
underlying forward rate. 

The key idea of the paper is to assume that the state-space of the realized variance process 
lies on a circle. This allows us to apply a block-diagonalization algorithm (described in detail 
in Appendix [X]) and the standard methods from spectral theory to obtain a joint probability 
distribution function for the underlying process and its realized variance or volatility at any time in 
the future (e.g. for time horizons of 6 months and 1 year see Figure [T5]) . This joint pdf is precisely 
what is needed to price a completely general payoff that depends on the realized variance and on 
the underlying. 

There are two natural and useful consequences of this approach. One is that we do not need to 
specify exogenously the process for the variance and then try to find an arbitrage-free dynamics for 
the underlying, but instead imply such a process from the observed vanilla market via the model for 
the underlying (a term-structure of the fair values of variance swaps as implied by the vanilla market 
data and our model is shown in Figure [TS]) . The second consequence is that this approach bypasses 
the use of Monte Carlo techniques and therefore yields sharp and easily computable sensitivities of 
the volatility derivatives to the market parameters. This is because the pricing algorithm yields, 
as a by-product, all the necessary information for finding the required hedge ratios. 

There is a rapidly growing interest in trading volatility derivatives in financial markets, mainly 
a consequence of the following two factors. On the one hand, pure volatility instruments are used 
to hedge the implicit vega exposure of the portfolios of market participants, thus eliminating the 
need to trade frequently in the vanilla options market. This in itself is advantageous because of the 
relatively large bid-offer spreads prevailing in that market. On the other hand, volatility derivatives 
are a useful tool for speculating on future volatility levels and for trading the difference between 
realized and implied volatility. 

This interest is reflected in the vast amount of literature devoted to volatility products. The 
analysis of the realized variance is intrinsically easier than that of realized volatility because of the 
additivity of the former. Under the hypothesis that the underlying price process is continuous the 
realized variance can be hedged perfectly by a European contract with the logarithmic payoff first 
studied in [26] and a dynamic trading strategy in the underlying. This approach does not require 
an explicit specification for the instantaneous volatility process of the underlying and can therefore 
be used within any stochastic volatility framework. This idea has been developed in [TT] and |14j . 
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where the static rephcation strategy for the log contract, using calls and puts, is described. Types of 
mark-to-market risk faced by a holder of a variance swap are studied and classified in [13] . A direct 
delta- hedging approach for the realized variance is given in [55] . A shortcoming of pricing variance 
swaps without specifying a volatility model (as described in [11] and [15]) is that this methodology 
does not yield a natural method for the computation of the sensitivities to market parameters (i.e. 
the Greeks). In [2^ a diffusion model for the volatility process is specified that allowed the authors 
to use PDE technology to price and hedge variance and volatility swaps, as well as more general 
payoffs. Derivatives on the realized volatility can be considered naturally as derivatives on the 
square root of the realized variance. In [B] the authors provide a volatility convexity correction which 
relates the two families of derivatives. Another approach, pioneered in [U], develops a robust hedging 
strategy for volatility derivatives which is analogous to the one for variance swaps. This method 
works for continuous processes only and is based on the observation that, under the continuity 
hypothesis, there is a simple algebraic relationship between the Laplace transform of the process and 
the Laplace transform of its quadratic variation. There are some technical difficulties in computing 
the relevant integrals, which have been dealt with in [18]. In [31] the authors investigate hedging 
techniques for discretely monitored volatility contracts that are independent of the instantaneous 
volatility dynamics. Another model-independent hedging approach in an environment with jumps 
for variance swaps is presented in [29]. A modelling framework, analogous to the HJM, for the 
term-structure of variance swaps has been proposed in [7]. The starting point is the specification of 
the function-valued process for the forward instantaneous variance, which yields an arbitrage-free 
dynamics for the underlying. In order to be calibrated, this model requires at time zero the entire 
variance swap curve. 

The paper is organized as follows. The key idea that allows us to price general derivatives on 
the realized variance is introduced in Section [2l Section [3] explains the pricing algorithm, based 
on Theorem 13.11 for derivatives on the realized variance. In Section H] we discuss the calibration 
of the model to a wide range of strikes and maturities for options written on the S&P 500. In 
Section [5] we carry out numerical experiments and consistency tests on the calibrated model, in- 
cluding a comparison with a Monte Carlo method (see Subsection 15. 6p . Concluding remarks are 
contained in Section [HI The numerical algorithm required to make this idea applicable is described 
in Appendix lAl which also contains the proof of Theorem 13.11 Appendix IB] describes some of the 
volatility contracts that can be priced within our framework and are referred to throughout the 
paper. 



2. The variance process for continuous-time Markov chains 

As mentioned in the introduction the main idea of this paper is to describe the realized variance 
of the underlying risky security as a continuous-time Markov chain on a state-space which lies on 
a circle. The model for the underlying Ft will be given in the forward measure as a mixture of 
local and stochastic volatility coupled with an infinite activity jump process. It will be defined 
on a continuous-time lattice in a largely non-parametric fashion, as described in [1], where it was 
used for modelling the foreign exchange rate. We will briefly recall the definition of the model in 
Section [3] and describe its calibration to the vanilla surface on the S&P 500 in Section [H 
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In the current section and the one that foHows our task is to build a mechanism that will make it 
possible to identify the random behaviour of the realized variance of the process Ff. The approach 
we take will only use the fact that the process Ft is a continuous-time Markov chain and will not 
depend on any other properties of the process. In order to do this we must return briefly to the 
fundamental theory (see [27], Chapters 2 and 3, for background on continuous-time Markov chains). 

Let L be a generator for a continuous-time Markov chain Ft defined on a finite state space Q. 
In other words the operator L is given by a square matrix (L{x, y))x,y£n which satisfies probability 
conservation and has non-negative elements off the diagonal. Each element L(2;, y) describes the 
first order change of the probability that the chain Ft jumps from level F{x) to the level F{y) in 
the time interval [t,t + dt], where the deterministic function F : — > (0, oo) is an injection which 
defines the image of the process Ft- 

Our aim in this section is to extend (i.e. lift) the generator L to the generator L, which will 
describe the dynamics of the lifting (Ft, It) of our original process Ft. The component It will be a 
finite-state continuous-time Markov chain, which will be described shortly, that approximates the 
realized variance, up to time t, of the underlying process Ft. The generator of the chain {Ft, It) 
will give us the probability kernel for the process It, which is what we are ultimately interested in 
as it contains the probability distribution of the relevant path information. 

We will start by describing a continuous process T,t which will contain the relevant path- 
information of the underlying Markov chain Ft and then define a finite-state stochastic process 
It that will be used to approximate Sj. The key feature of the state-space of the process It, which 
allows numerical tractability of our approach, is that it is contained on a circle in the complex 
plane (see Section 

The procedure we are about to describe works for a specific type of path-dependence only. 
Assume that if at time t the underlying process Ft is at a state F{x) for some x £ ft, then the 
change of the value St, in the infinitesimal time interval dt is of the form 



where the function Q : Q ^ M, defined in terms of the underlying stochastic process, has two key 
properties: 

• the mapping Q is independent of the path the underlying process follows on the interval 
[0,t) and 

• the value Q{x) depends on the level x at time t and on the distribution of the underlying 
process in the infinitesimal future time interval dt as given by the Markov generator L. 

The first property states that the change in the path- information St, over a short time period 
[t,t + dt], does not depend on the path taken by the underlying process up to time t. The second 
property tells us that the evolution of the path information over the interval [t,t + dt], conditional 
on the current level of the process Ft, is determined by the future distribution of Ft over the 
infinitesimal time interval. 

It is clear that the realized variance of the process Ft up to time t, can be captured by a random 
process Sj. Indeed, if we define the function Q in the following way 



(iSt = Q{x)dt, 



(1) 
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it follows immediately that the above conditions are satisfied. This is because the realized variance 
of Ft is simply an integral over time of the instantaneous variance of Ft which is given by ([1]) . 

The last observation is crucial for all that follows, because it implies that the process St is 
uniquely determined by its state-dependent instantaneous drift. In particular we see that the pro- 
cess Tit has no volatility term and that it has continuous sample paths since it allows a representation 
as an integral over time. 

The fundamental consequence of these facts is the following: a finite-state Markov chain It can 
be used as a model for the process if and only if the instantaneous drift of the chain is equal to 
Q{x), whenever the underlying process Ft is in state F{x), and the instantaneous variance of It is 
equal to zero up to first order. The first requirement clearly follows from the discussion above. The 
second condition is there to reflect the fact that the process Sj is a continuous Ito process which 
has no volatility term. A non-constant random process on a lattice will always have a non-zero 
instantaneous variance, but the second condition ensures that this instantaneous variance goes to 
zero as quickly as the lattice spacing itself. In other words the Markov chain approximation It for 
the process Sj must exhibit neither diff'usion nor jump behaviour. Therefore, as we shall see in the 
next subsection. It must be a Poisson process with state-dependent intensity. 

2.1. Lifting of the generator for the underlying process. Let It as above denote the Markov 
chain whose value approximates the realized variance of the forward price Ft from time to 
time t. We shall express It as arrit where nit is a Poisson process (with non-constant intensity) 
starting at and gradually jumping up along the grid given by = {0, . . . ,2C}, where C is an 
element in N. The constant a controls the spacing of the grid for the realized variance It- 

We are now going to specify precisely the dynamics of the process nit which is a fundamental 
ingredient of our model. As mentioned before, the process nit will behave as a Poisson process 
whose intensity is determined by the level of the underlying. In other words the Markov generator, 
conditional on the underlying process being at the level F{x), is of the form 



' ^Q{x) if d = (c + 1) mod (2C + 1); 

(2) L'^ix : c,d) := i 

, -iQ(x) if d = c. 

The variables c and d are elements of the discrete set ^ = {0, ...,2C} and the function Q is 
the instantaneous variance of the underlying process as defined in ([T]). This family of generators 
specifies the dynamics of the process It with the following instantaneous drift 

2C 



lim E 

Ai-»0 



h+At — h 



At 



Ft = F{x),It = I{c) 



^{I{d)-I{c))L^{x:c,d) 



d=0 



(3) 



1 



a{c + 1 — c)—Q{x) = Q{x 



a 



for all values F{x) of the underlying process Ft and all integers c which are strictly smaller than 
2C. This implies that the Markov chain It has the same instantaneous drift as the actual realized 
variance. A similar calculation tells us that the instantaneous variance of It equals aQ{x). Since 
a is the spacing of the lattice for the chain It, we have just shown that the first two instantaneous 
moments of It match the first two moments of the realized variance for all points on the lattice 
{0, . . . , 2C} except the last one. 
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Notice however that the equahty ([3]) breaks down if c equals 2C, because we have imposed 
periodic boundary condition for the process It- Put differently this means that It is in fact a 
process on a circle rather than on an interval. The latter would be achieved if we had imposed 
absorbing boundary conditions at 2C. That would perhaps be a more natural thing to do since the 
process Sf is certainly not periodic. But an absorbing boundary condition would destroy the delicate 
structure of the spectrum of the lifted generator L which is preserved by the periodic boundary 
condition. It is precisely this structure that makes the periodic nature of It a key ingredient of our 
model, because it allows us to linearize the complexity of the pricing algorithm (see appendix It 
should be noted that the general philosophy behind either choice of the boundary condition would 
be the same: the lattice in the calibrated model should be set up in such a way that the process 
never reaches the boundary value 2C, because if it does an inevitable loss of information will ensue 
regardless of the boundary conditions we choose. 

We are now finally in a position to define the lifted Markov generator of the process (Ft, It) as 

(4) L(x, c; y, d) := L{x, y)6c4 + ^"^{x : c, d)6x,y, 

where x,y are in the state-space of the chain Ft and c,d are elements in ^. The structure of the 
spectrum of the operator L will be exploited in Section [3] to obtain a pricing algorithm for payoffs 
which are general functions of the realized variance (see Theorem 13. ip . The reason for specifying 
the generator Ij"^(x : c,d) by ([2|) (i.e. insisting on the periodic boundary condition for the process 
It) will become clear in the proof of Theorem 13.11 which exploits the spectral properties of the 
lifted generator L (see appendix |X] for the precise description of the spectrum of L). 

3. Pricing of derivatives on realized variance 

Recall that the generator L for the underlying continuous-time Markov chain Ft in the forward 
measure, described in [I], is defined on the state-space QxV via an injective function F : 0,xV ^ 
(0, oo). The set il. with elements is the state-space of the various jump-diffusions and the set 
V with M elements is the state-space for the random switching between regimes (see the formulae 
on page 7 in [1] for the precise definition of the generator matrix (L(a;, /3; y, 7))(x,/3),(y,7)gr2xy )■ 
Section S] we shall see that for the options data on the S&P 500 used in this paper, it suffices to 
take N = 76 with only two regimes (M = 2). 

Let TiT be the realized variance over the time interval [0, T], expressed in annual terms, as defined 
in (I16p . In this section we are going to find pricing formulae for general payoffs that depend on the 
annualized realized variance Yjj^. 

Our first task is to find the probability kernel of the lifted process (Ft, It) which was defined in 
Section O Recall that the Markov chain It = arrit , which is used to model the realized variance, 
is specified in terms of a translation invariant process mt, whose domain is ^' = {0, . . . , 2(7}, and 
some positive constant a which determines the lattice spacing for the domain of It- The value 
C G N specifies the size of the lattice for the realized variance and has to be large enough so that 
the process It, starting from zero, does not transverse the entire lattice. This is a very important 
technical requirement as it ensures that there is no probability leakage in the model (which is 
theoretically possible since we are using periodic boundary conditions for the process It). The 
dynamics of the chain nit are given by the Markov generator in ^ . 
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Recall also that the process It records the total realized variance up to time t, which implies 
that the annualized realized variance that interests us, will be described by ^It where the time 
horizon T is expressed in years. The key ingredient in the calculation of the probability distribution 
function of the process (Ft, It) is the block-diagonalization algorithm from subsection IA.4I of the 
appendix. We are now going to apply it to the generator ^ in order to find the joint pdf of the 
lifted process. 

3.1. Probability kernel of the lift (Ft, It). The generator L(x, /3, c; 7, d) of the process (F^, It), 
given by ([1]) , is a square partial-circulant matrix as defined in appendix El It is also clear that 
the dimension of the vector space acted on by the matrix L is MN{2C + 1). The coordinates 
(x,/?, c), {y,j,d) of the matrix L (i.e. the lattice points of the process) lie in the set Q x V x ^ , 
where is the grid for the underlying forward rate and the set V contains all volatility regimes 
of the model for the forward Ft. Notice that the circulant matrices L™ from ([2|), used in the 
definition of L, are very simple because the only non-zero elements are on the diagonal, just above 
the diagonal and the element in the bottom left corner of the matrix. In other words if we interpret 
the matrix L™, associated with the lattice point {x,(3), in terms of the definition of a circulant 
matrix given at the beginning of appendix [Aj we see that 

Q{x,P) 

Cl = -Co = , 

a 

where the function Q{x,(3) is the instantaneous variance as defined by ([T|) and the constant a is 
the lattice spacing of the domain of It- All other elements Cj are equal to zero. 

It is therefore clear that, using the expression ([8]) for the eigenvalues of circulant matrices, 
equation (|15p can be reinterpreted as 

Q{x,j3) 

Lfc(2;, /?; y, 7) := H^, /?; 7) + \x,i3),{yn) (^~*^'° ~ 1) — — ' 

where is the A;-th block in the block-diagonal decomposition of L and the value of is given by 
the expression 

(5, 

The index k in these expressions runs from to 2C. Notice that the matrices L^. differ from the 
Markov generator L only along the diagonal. We are now in a position to state the key theorem 
that will allow us to find the probability kernel of the lifted process. 

Theorem 3.1. Let L be the Markov generator of the stochastic process {Ft, It) as described in 
Section and let (p : C ^ C be a holomorphic function. Then the following equality holds 

1 

</.(i)(x,/3,c;y,7,d) = ^^—Y,e-''"^""'U{Lk){x,P;y,j), 

k=0 

where Lk is the operator defined above, pk is given by ^ and (x, /?, c), (y, 7, d) are elements of 
QxV x^. 

Before embarking on the proof of this theorem we may summarize as follows: if a linear operator 
A can be block-diagonalized by a discrete Fourier transform (cf. last paragraph of subsection IA.4I 
in the appendix), then so can any operator 4>{A) where i;^> is a holomorphic function defined on 
the entire complex plane. Note that the assumptions on function (j) and the linear operator L in 
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Theorem 13.11 are too stringent and in fact the theorem holds in much greater generaUty. For our 
purposes however the setting described in the theorem is sufficient as it appUes directly to the 
model. Therefore the proof of the theorem, given in appendix IA.51 only applies in the restricted 
case stated above. Note that the idea of semi-analytic diagonalization, which is central to the proof 
of our theorem, has been applied extensively in physics and numerical analysis (see for example [21] . 
[2], [T7] and the references therein). 

We can now find the full probability kernel of the process (Ft, It) by applying Theorem 13.11 
in the following way. More explicitly, for any pair of calendar times t and T, such that t < 
T, let the conditional probability F{Ft = F{y,^),lT = d\Ft = F{x,(3),It = c) be denoted by 
p({x, (5, c),t; {y, 7, d),T). The well-know exponential solution to the backward Kolmogorov equation 
and Theorem 13.11 imply 

= J—i E e-^^(-'^)e(/(^)-/W)I^^ (X, P- y, 7) 

k=0 

2C M-N 

(6) = ^5^J]e^^(/(^)-/W)e--(-'^)<(x,/3).^(,,7), 

fc=0 n=l 

where are the eigenvalues and v!^ are the eigenvectors of L^. As usual we denote by the 
columns of the matrix U'j^^ , where Uk consists of all the eigenvectors of L^. Function / in the above 
formula is the deterministic time-change for the underlying model which was introduced in [1], 
subsection 2.4 (see Figure [3] for the specification of function / implied by the volatility surface for 
the S&P 500 used in Section g]). 

Formula ([6]) is our key result because it allows us to price any derivative which depends jointly 
on the realized variance of the underlying index and the index itself at any time horizon T. In the 
following subsection we state an explicit formula for the value of such a derivative. 

3.2. Pricing derivatives on the realized variance. Let us assume that we are given a general 
payoff h(TiT-t) that depends on the annualized realized variance St-* between the current time t 
and some future expiry date T. The current price of a derivative with this payoff within our model 
can be computed directly using the joint probability distribution function ^ in the following way 

Ct{x, P) = e-('-(T)T-r-(t)t) ^ I ^ ^((^^ 0), t; {y, ^,d),T)] h 

The sum in the brackets is the marginal distribution of the process It at time T — t. Notice that 
the realized variance process must always start from at the inception of the contract. The factor 
normalizes the value of Ix-t so that it is expressed in annual terms. As before, the constant a 
is the lattice spacing for the realized variance. 

The point x from is chosen in such a way that the equation F{x) = e^^^^^~'^^^^^^ St holds, where 
St is the index level at the current time t, the functions r, d : [0, T] — > M4. are the deterministic 
instantaneous interest rate and dividend yield respectively and /? in V corresponds to the volatility 
regime at time t. 

Notice that because in the valuation formula above there is no restriction on the payoff function 
h, the contract Ct can be anything from a variance swaption to a volatility swap (or swaption) 



f ad \ 

Vr^t 
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and can be priced very efficiently using the calibrated model. Since formula ^ yields the joint 
distribution of the realized variance '^T-t and the underlying security St, a trivial modification 
of the above expression would give us the current price in our model of any payoff of the form 

h{ST, ^T-t)- 

4. Calibration to the vanilla surface 

We are now going to calibrate the model for the underlying, described in Section 2 of p], to 
the implied volatility surface of the S&P 500 equity index for maturities between 1 month and 5 
years (see legend of Figure [T] for all market defined maturities used in the calibration) and a broad 
range of liquid strikes for each maturity. The market data consist of the implied Black-Scholes 
volatilities for each strike and maturityO Our task is to find a set of values for the parameters of 
the model such that if we reprice the above options and express their values in terms of the implied 
Black-Scholes volatilities, we reobtain the market quotes. 

The first choice we need to make pertains to the number of regimes M we are going to use. In 
foreign exchange markets the implied volatility skew exhibits complex patterns of behaviour across 
currency pairs and the corresponding model required no fewer than 4 regimes to accommodate this 
complexity (see [T]). In equity markets, by contrast, the skew is always large and negative and 
becomes less pronounced as time to maturity increases (see chapter 8 in ^9\). Thus in the case of 
S&P 500 two regimes prove sufficient (M = 2). The number of regimes cannot be smaller than two 
because in the case of a single regime the model becomes a jump diffusion which, as is well-known, 
does not describe adequately the implied volatility surface (see chapter 5 in [T9]). 

In order to calibrate our model we select an non-homogeneous grid with 76 = points used to 
span the possible values of the forward rate Fp. The grid is elliptical in the following sense: close 
to the current value of the spot, the lattice points are very dense; the density goes down as we go 
further away from the current level of spot. The top node is at 100^ and the bottom node at 100~^. 

The underlying grid for the model is therefore of the size MN = 2 • 76 = 152, which makes the 
model more efficient computationally than the one used for modelling the foreign exchange rate, 
where the lattice was roughly double in size. 

The next group of model parameters which do not need to be changed with every calibration are 
the levels Fq and Fi and the stochastic volatility matrices and QY • Since the model is defined 
so that the starting regime is regime 1, we take the level Fi to be equal to the current value of spot 
100. The level Fq is chosen somewhat arbitrarily to be 95. If the underlying starts trading around 
this level, regime assumes a more dominant role in the model. If we took Fq to be, for example, 
90, that would have a minor effect on the short term (up to six months) implied volatility skew of 
the model and a negligible effect on longer maturity skews. This is because by moving Fq from 95 
to 90 we decrease the influence of regime 1 on the short term maturity skew. 

The stochastic volatility generators Qq and QY are chosen to reflect the fact that in equity 
markets, a downward move in the underlying index results in steeper negative skews. The starting 
regime 1 corresponds to the current skews for short maturities, whereas regime corresponds to 
steeper negative skews, which will occur if the value of the underlying index drops. The generator 



Notice that we are using more strikes for longer maturities than for shorter maturities. This is not an inherent 
requirement but a consequence of the initial structure of the market data we used for calibration of the model. 
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Qq , which is dominant if the index is trading around 95 or below, assigns a very shght probabihty 
to changing back to regime 1. This is because if the index is trading at 95 or below, the skew is 
unlikely to become less steep. For the same reason assigns a substantial probability to staying 
in regime 0. On the contrary, the generator QY assigns roughly the same probability to staying 
in regime 1 and to switching to regime 0. This is because when the underlying index is trading 
close to the current level of spot, the skews can become less or more steep with roughly the same 
probability. 

In order to calibrate to the specific instance of the volatility surface, the parameters in the CEV 
jump diffusion need to be chosen. Note that we can calibrate the model to the entire implied 
volatility surface without recourse to the parameters CTq, for a = 0, 1. This is due to the fact that 
in our model they only affect option values for strikes below 20. Since, when calibrating to the 
implied volatility surface, we are not interested in strikes that are so far away from the at-the- 
money strike, we could take the value of the parameters cTq, for a = 0, 1, to be 10000% (this choice 
amounts to excluding the paramters cJq, for a = 0, 1, from the model). 

On the other hand, the values of variance swaps are very sensitive to the strikes below 20 (see 
expression (|20|) in appendix lB.3l for the hedging portfolio consisting of vanilla options) and therefore 
exhibit a strong dependence on the CTq,. Therefore if we set cJq, = 10000%, for a = 0, 1, variance 
swaps would be grossly overpriced. Since parameters CTq,, for a = 0, 1, do not affect the at-the- 
money skew, we could use these parameters to calibrate the model to the term structure of variance 
swaps. In this paper we did not have the relevant variance swap data to calibrate to, we therefore 
chose ctq = ai = 60% which gave reasonable values for the term structure of variance swaps (see 
Figure I18|) . Note also that ai influences more markedly the short term variance swaps and cjo has 
a greater effect on longer term contracts. Therefore the steepness of the term structure of variance 
swaps implied by the model can be adjusted by taking different values for cji and o-q. 

Going back to the remaining parameters of the CEV jump diffusion, note that, since there are 
no positive jumps in the equity index, we take f^, for a = 0, 1, to be 0. We take the CEV volatility 
o"! in regime 1 to be close to the at-the-money implied volatility for the shortest maturity we are 
calibrating to (in this case one month). The parameters /3i and are chosen to fit the skew for 
the maturity of one month. 

Introducing a second regime into the model does not alter the one-month skew (this follows 
from the definition of the stochastic volatility generator GY)- Since regime corresponds to the 
underlying index trading at a lower level, the values of /9o and i^q are more extreme, as they 
correspond to the tightening of the skew. Their values are chosen in such a way that the model 
qualitatively fits the skews for all maturities (i.e. the risk-reversals are priced correctly but the 
level of the at-the-money implied volatility is not necessarily matched). 

By choosing the parameters using these guidelines (see tabled]), we obtain an implied volatility 
surface that has the correct qualitative shape, but the at-the-money levels are not represented 
correctly. To remedy this, mild explicit time-dependence needs to be introduced. This is achieved by 
introducing a deterministic function which transforms calendar time to financial time (see Figure [3]). 
The effect on the implied volatility skew is that of a parallel shift for each maturity. In other words, 
the shape of the skew remains unchanged. For the final result see Figure [TJ 
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The main calibration criterion has been to minimize the exphcit time-dependence in order to 
preserve the correct (i.e. market impHed) skew through time. The stationarity requirement in the 
cahbration ensures that the forward ske^v^l^ for any maturity retains the desired shape as can be 
observed in Figures El El [9] and [lOl 

The S&P 500 index equaled 1195 at the time when the option data were recorded. Throughout 
the paper a relative value of the index, set at 100, is used for simplicity. The forward price levels 
F(x), where x is an element of the underlying grid $7, are also measured on the relative scale. As 
mentioned above our starting volatility regime is regime a = 1. 



a 




f3a 


aa 






Fa 





16% 


-0.8 


60% 


0.18 





95.00 


1 


13% 


-0.3 


60% 


0.15 





100.00 



Table 1. Parameters for the local volatility regimes and the jump intensities. 



Markov generators for stochastic volatility (a = 0, 1). 

The code used to obtain these results is written in VB.NET and relies on LAPACK for the 
computation of the spectrum of the generator. In the current model the generator is a matrix of 
the size 152 x 152. The time required to diagonalize the generator and price 76 options per maturity 
for 6 maturities (i.e. 456 options) is less than 4 seconds on a single Pentium M processor with 1GB 
of RAM. Note that the time consuming step is the computation of the spectrum and eigenvectors. 
Once these are known the probability kernel for any maturity, and therefore all option prices for 
that maturity, can be obtained within a fraction of a second. 

4.1. The parameters C and a in the lift {Ft, It)- Figures [T5 l 1161 and 1 171 contain the graphs of 
the joint distribution functions of the spot level and the annualized realized variance for maturities 
between 6 months and 5 years. The marginal distributions of the annualized realized variance, 
obtained from the above joint distributions by integration in the dimension of the spot value of the 
index, are shown in Figure [TTl 

The parameters in the model that influence the dynamics of the realized variance It are the 
number of lattice points C and the lattice spacing a. It turns out that for maturities up to 5 years 
it is sufficient to take 201 = 2C + 1 uniformly distributed points for the grid of the realized variance 
(i.e. C = 100). Table [2] contains the values of a for different maturities. It should be noted that 
the numerical complexity of the algorithm used to obtain the pdf of the joint process (Ft, It) for 
different maturities is constant since C does not change. As mentioned earlier, it is crucial that 
the value of C be chosen large enough with respect to the spacing a so that the process cannot get 

^There is a closed form solution for the value of a forward-start in the Black-Scholes model and the only unknown 
parameter in that formula is the forward volatility. It is therefore customary to define a forward smile of any model as 
a function mapping the forward strike to the implied forward volatility which is obtained by inverting the Black-Sholes 
formula (see subsection IB. 31 for the precise definition of a forward strike) . 
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6m 




2v 


3v 


4v 


5v 


a 


0.00085 


0.00165 


0.00362 


0.00570 


0.00759 


0.00970 



Table 2. Parameters specifying the geometry of the lattice for the reahzed variance 
process It- 

to the other side of the lattice with positive probability. The guiding principle has been to ensure 
that the probability of this event is smaller than 10~^. In fact for a fixed maturity t up to five 
years the formula 

49^ 

gives a definition of the lattice spacing a{t) which has the required property (for the graph of a(t) 
see Figure [Til) . Here the function f{t) is the financial time introduced in [1] (see Figure [3] for the 
graph of /(t) used in calibration to the volatility surface). The intuition behind this definition 
is that the maximal value of the realized variance increases linearly with financial time and is 
independent of the number of lattice points used to model it. With this choice of parameters the 
process It wraps around with probability lO"*^ if t equals 6 months, and with probability 10~^ 
when i is 5 years (these probabilities can be made even smaller -at the expense of computational 
efficiency- if we choose larger values for C and smaller spacing a(t)). 

The computational time required to obtain the joint probability law of {Ft, It) for a fixed time 
t on a Pentium M processor with 1GB of RAM, using VB.NET and relying on LAPACK for the 
calculation of the spectra, is no more than 130 seconds. The task would be to diagonalize 201 
complex matrices of the size 152 x 152. However the symmetry of the discrete Fourier transform 
allows us to reduce the number of matrices we need to diagonalize to 101. Since these matrices 
are independent of each other, the algorithm for obtaining the joint law is parallelizable and the 
computational time can be reduced significantly by using multithreaded code running on several 
processors. 

5. Numerical results 

Let us now use the calibrated model together with the pricing and hedging algorithms described 
in [l], subsection l3.2l and appendixElto perform some numerical experiments and consistency checks 
for vanilla options, forward-starts, variance swaps and other volatility derivatives. We conclude this 
section by a numerical comparison of the spectral method developed in Sections [2] and [3] with a 
Monte Carlo pricing algorithm for volatility derivatives (see subsection 15. 6p . 

5.1. Profiles of the Greeks. In order to test the pricing methodology of the model for the 
underlying forward rate we pick a strip of call options with the same notional but with varying 
maturities, all struck at the current spot level. Since the entire framework for the underlying is 
expressed in relative terms with respect to the current value of the index, the strike used for this 
strip of options is 100. 
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Because we are interested in the behaviour of our pricing algorithm in changing market condi- 
tions, we are going to study the properties of delta, gamma and vega, given by the formulae 

Co(x + l,l)-C7o(x-l,l) 



A(x) 

r(x) 



F{x + 1) - F{x - 1) ' 
^ Co(x + 1, 1) + Co(x - 1, 1) - 2Cq{x, 1) 

{F{x + l)- F{x-l)Y 
Co(x,0)-C7o(x,l) 



(To - CFi 

respectively, as functions of the current level of spot. Here x denotes the point in that corresponds 
to the spot level of the index at time 0. 

This task does not pose any additional numerical difficulties because all it requires is the knowl- 
edge of the probability distribution functions of the underlying at the relevant maturities conditional 
upon the starting level, which can be any of the points in the grid But if we have already priced 
a single option, then these pdfs are available to us without any further numerical efforts. This is 
because, when pricing an option, the algorithm described in [1], Section 3, calculates the entire 
probability kernel for each starting point in the grid, even though the option pricing formula only 
requires one row of the final result. In this situation we require all the rows of the probability ker- 
nel so as to obtain the prices of our option, conditional upon different spot levels, by applying the 
matrix of the kernel to the vector whose coordinates are values of the payoff calculated at all lattice 
points. As defined by the formulae above, the Greeks are linear combinations of the coordinates of 
the final result of the calculation. 

Figure H] contains the delta profile of the call options for different maturities. Figures [S] and 
give the gamma and the vega profiles of the call options in the strip. As expected, an owner of a 
vanilla option is both long gamma and long vega. The shapes of the graphs in Figures [5] and [6] also 
confirm that, according to our model, the at-the-money options have the largest possible vega and 
gamma for any given maturity. A cursory inspection of the scales of gamma and vega for options 
with the same notional indicates that in our model some calendar can be simultaneously 

long vega and short gamma or vice versa. 

5.2. Forward smile. The forward volatility can be defined using the Black-Scholes formula for 
the forward-starting options (see appendix IB. 31 for details). One of the parameters in this formula 
is a forward strike. Hence any model for the underlying process defines a functional relationship 
between forward strikes and implied forward volatilities using the Black-Scholes pricing formula 
for forward-starts, in much the same way as it defines the implied volatility as a function of strike. 
This functional relationship is known as the forward smile. 

The reason why the forward smile is so important lies in the fact that it determines the conditional 
behaviour of the process. It is well-known that knowing all the vanilla prices (i.e. the entire implied 
volatility surface) is not enough to price path-dependent exotic options. In terms of stochastic 
processes this statement can be expressed by saying that knowing the one-dimensional marginals 
of the underlying for all maturities does not determine the process uniquely. The forward smile 
contains the information about the two-dimensional distributions of the underlying process. In 



''A calendar spread is a structure defined by going long one call option and going short another call option of the 
same strike but different maturity. 
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other words one can have two different models that are perfectly calibrated to the implied volatility 
surface but which assign completely different values to the forward-starting options (see e.g. |30]). 

Market participants can express their views on the two-dimensional distributions of the under- 
lying process by setting the prices of the forward-starting options accordingly. It is therefore of 
utmost importance for any model used for pricing path-dependent derivatives to have the implied 
forward smiles close to the ones expected by the market. Since the market implied forward volatil- 
ities are rarely liquid, the main market expectation is that the future should be similar to the 
present. Figures O [H [9] and [10] contain the forward smiles implied by our model for maturities 
between 3 months and 2 years. Since we have no market data to compare them with, we can only 
say that the qualitative nature of the implied forward smiles is as expected in the following sense: 
for a fixed time T' the forward smiles are flattening with increasing T (for definition of times T' 
and T see appendix IB.Sh and for a fixed difference T — T' the shapes of the forward smiles look 
similar when compared across maturities T' . It should be noted that the latter point exemplifies 
the stationary nature of the underlying model and shows that we did not have to use extreme values 
of the model parameters to calibrate it to the entire implied volatility surface, because the two- 
dimensional distributions of the underlying process have the stationary features which are expected 
by the market participants. 

We should also note that a statistical comparison of the forward volatility smiles implied by 
the model and by the market can be carried out easily because pricing a forward-starting option 
consists of two consecutive matrix-matrix multiplications which require little computing time. 

5.3. Probability distribution for the implied volatility index. One of our goals in this paper 
has been to describe the random evolution of VIX through time. In Figured^we plot the probability 
density functions of the volatility index, as defined at the end of appendix IB. 31 

Since we can calculate explicitly distributions of VIX for any maturity (like those in Figure [T9]) , 
pricing a European option on the VIX in our framework amounts to summing the values of the 
payoff against the pdf in the same way as was done for European payoffs on the underlying security 
in [I]. 

Note also that the calculation of the distributions of VIX is independent of the lifting procedure 
described in Sections and O Therefore the computational effort required to obtain it is minimal 
as the task at hand solely consists of pricing portfolios of forward-starting options. 

5.4. Distribution of the realized variance. Most of the modelling exercise presented so far has 
been directed towards finding the probability distribution function for the annualized realized vari- 
ance of the underlying process. Figure [TT] contains the pdfs for the realized variance for maturities 
between 6 months and 5 years as implied by the model that was calibrated to the market data 
for the S&P 500. The corresponding term structure of the fair values of variance swaps for these 
maturities is given in Figure [TSl This term structure is obtained by calculating the expectation (in 
the risk-neutral measure) of the probability distribution functions of Figure [TTl 

Figure [18] also contains the current value of the logarithmic payoff, as given by (I20p . for the above 
maturities. In [15] the performance of the log payoff as a hedge instrument for the variance swap 
is studied in the presence of a single down-jump throughout the life of the contract. In equation 
42 the authors show that, in this case, the log payoff is worth more than the variance swap. Note 
that the values of the log payoff in Figure [TH] dominate the values of the variance swaps, as given 
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by our model. This agrees with [15] because our model for the underlying exhibits jumps down 
but no up-jumps. In chapter 11 of [19] it is shown that for a compound Poisson process with 
intensity A = 0.61 and normal jumps with mean —0.09 and variance 0.14, the difference between 
the current value of the logarithmic payoff and the variance swap (both maturing in one year) 
equals 0.3% if both prices are expressed in volatility terms (i.e. the difference between the square 
roots of the values of the annualized derivatives). This should be compared to 0.5355% which is 
the corresponding difference in our model (see table [3]). The prices, expressed in volatility terms, 
seen in Figure [18] are given by the following table: 



Maturity T 




l°g(t) 




Portfolio (fTHD 


Eo [Et] 


Eo 






0.5 


15.4613 


15.4620 


15.0054 


14.6960 


1 


16.0833 


16.0836 


15.5478 


15.0378 


2 


18.6775 


18.6776 


17.9255 


16.7808 


3 


20.8223 


20.8224 


19.8931 


18.1161 


4 


22.0701 


22.0703 


21.0327 


18.8013 


5 


23.5183 


23.5186 


22.3735 


19.6923 



Table 3. The random variable is the annualized realized variance of the under- 
lying process. All the prices are quoted in volatility terms. Figure [TSl contains the 
pictorial description of the data presented in this table. 

The class of HJM-like models for volatility derivatives (see for example [7]) require an entire 
term structure of variance swaps (like the one in Figure [TS|) to be calibrated. Our approach on the 
other hand implies one from the vanilla options data (and of course some modelling hypothesis). 
However, as mentioned in Sectiond] the parameters cJo, for a = 0, 1, have a negligible effect on the 
model implied volatility of vanilla strikes that are observed in the markets but have a substantial 
effect on the option prices for strikes smaller than 20. It is clear that portfolio (I18p . and therefore 
the corresponding variance swap itself, is very sensitive to these options. Therefore by adjusting 
a a, for a = 0, 1, we can calibrate to variance swaps for at least two maturities without distorting 
the implied volatility surface of our model for strikes that are closer to the at-the-money strike. 
In particular the steepness of the variance swap curve implied by the model can be adjusted by 
choosing different values for cji, which controls the short maturities since regime 1 is the starting 
regime, and ctq, which has more of an influence on the far end of the curve. 

5.5. The log contract and variance swaps in a market without jumps. It is well-known 
that in a market where the underlying follows a continuous stochastic process the fair value of a 
variance swap is equal to the value of the replicating European option with the logarithmic payoff 
given in (f^U]) of appendix IB. 31 (see for example [K] or [H]). In the presence of jumps this equality 
ceases to hold as can be observed in Figure [TSl Recall from Section H] that in order to calibrate our 
model we had to use a non-zero value for the intensities of the down-jumps. 

In this subsection the aim is to confirm that the prices of variance swaps and logarithmic payoffs 
agree in our framework if the randomness of the underlying model is based purely on stochastic 
and local volatility. To that end we set up a simplified version of the model, with two volatility 
regimes only, using the following parameters: 
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^ OL 











10.0% 


0.70 


60% 








100 


1 


13.5% 


0.50 


60% 








110 



Table 4. Parameters for the local volatility regimes in the simplified model without jumps. 

/-CS 0.5 \ 
\ 0.5 -0.5/ 

Markov generators for stochastic volatility (a = 0, 1). 

Note that the jump intensities in this version of the model are deliberately set to zero. The 
corresponding probability distribution functions for maturities between 6 months and 5 years, as 
implied by the simplified model, are given in Figure [T2l 

We priced variance swaps, volatility swaps, logarithmic payoffs and structured logarithmic payoffs 
given by a portfolio of vanilla options in ()18l) of subsection IB. 31 on the notional of one dollar for 
maturities between 6 months and 5 years. The results, expressed in volatility terms, can be found 
in table [3 



Maturity T 




l°g(f) 




Portfolio (fT8]) 


Eo [St] 


Eo 






0.5 


10.3691 


10.3753 


10.3699 


10.3031 


1 


10.6368 


10.6450 


10.6387 


10.5327 


2 


10.9882 


11.0003 


10.9890 


10.8832 


3 


11.2036 


11.2142 


11.2013 


11.1186 


4 


11.3494 


11.3614 


11.3435 


11.2602 


5 


11.4574 


11.4691 


11.4513 


11.3593 



Table 5. The random variable Sj- is the annualized realized variance of the under- 
lying process. All prices are quoted in volatility terms. 

It follows from table [S] that the difference between the value of the logarithmic payoff and the fair 
value of the variance swap, according to our model, is less than 1 volatility point for all maturities. 
We can also observe the quality of the approximation of the portfolio of options (1180 to the log 
payoff as well as the convexity effect for volatility swaps. 

If the CEV parameters [3a, for a = 0, 1, take more extreme values, as they do in the calibrated 
model, the discrepancy between the current value of the log payoff and the fair value of the variance 
swap increases slightly even if jumps have not been introduced explicitly. We illustrate this by 
choosing the same model parameters as in table [T] of Section H] with the exception of i/qT, for 
a = 0, 1, which are taken to be 0. Table [6] shows the values of the relevant securities obtained using 
this model. 

The difference in value between the logarithmic payoff and the variance swap, both maturing in 
5 years, is less than 0.05%. This should be compared with the difference 0.0061% which occurs if 
the CEV parameters are less extreme (see table [5]). This phenomenon is due to the fact that any 
random move on a lattice is by its definition a (perhaps small) jump. If local volatility function is 
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Maturity T 


-|Eo 


l°g(t) 




Portfolio (fT8]) 


Eo [St] 


Eo 


v Sy 




0.5 




14.9461 


14.9437 


1 A fi77n 


1 




15.4234 


15.4181 




2 


1 7 snn9 

i / .oUUZ 


17.8003 


17.7871 


10. 1 / OU 


3 


1 Q 871 7 


19.8717 


19.8510 


1 8 1 791 


4 


21.1100 


21.1101 


21.0827 


18.9246 


5 


22.5537 


22.5538 


22.5171 


19.8770 



Table 6. The random variable T^t is the annuahzed realized variance of the un- 
derlying process and the full model with no jumps (i.e. u~ = 0, for q = 0, 1). All 
prices are quoted in volatility terms. 

steeper, such moves occur more frequently thus contributing to the larger difference. Note however 
that the density functions of the realized variance across maturities of the last example without 
jumps (see Figure [T3]) are qualitatively very close to the densities of the realized variance in the 
calibrated model (see Figure [TT]) . 

5.6. Comparison with Monte Carlo. In this subsection we perform numerical comparisons 
between the spectral method for pricing volatility derivatives, developed in Sections [2] and El and a 
direct Monte Carlo method. We generate 10^ paths for the underlying with a maturity of 5 years 
using the calibrated model (for parameter values see Section [3|), thus obtaining an approximation for 
the distribution of the realized variance for a set of maturities (1, 3 and 5 years). Each path consists 
of one point per day. Using this distribution we price the following derivatives: a variance swap 
■y/Eo [St], quoted in volatility terms, a volatility swap Eq [-v/Sy] and an option on the realized 
variance Eq [maxjST — (aii'o)^) 0}] , where the strike Kq := y^Eq [Sy] equals the current risk- 
neutral mean of the realized variance and the paprameter a equals 80%, 100%, 120%. The results 
are given in the following three tables. 



Maturity T = 5y 


Time (seconds) 


^Eo [St] 


Eo y^] 


a = 80% 


a = 100% 


a = 120% 


MC: 5 • 10^ paths 


402 


22.48% 


19.20% 


2.83% 


2.32% 


1.89% 


MC: lO^paths 


780 


22.46% 


19.19% 


2.81% 


2.31% 


1.88% 


Spectral method 


130 


22.37% 


19.69% 


2.61% 


2.05% 


1.59% 



Maturity T = 3y 


Time (seconds) 


VEo [St] 


Eo y^] 


a = 80% 


a = 100% 


a = 120% 


MC: 5 • 10^ paths 


219 


19.95% 


17.51% 


2.05% 


1.61% 


1.28% 


MC: 10^ paths 


439 


19.92% 


17.49% 


2.04% 


1.60% 


1.27% 


Spectral method 


130 


19.89% 


18.12% 


1.85% 


1.34% 


0.99% 



Maturity T = ly 


Time (seconds) 


^Eo [St] 


Eo 


\/St 




a = 80% 


a = 100% 


a = 120% 


MC: 5 • 10^ paths 


65 


15.53% 


14.16% 


1.13% 


0.80% 


0.56% 


MC: 10^ paths 


130 


15.53% 


14.15% 


1.13% 


0.80% 


0.56% 


Spectral method 


130 


15.55% 


15.04% 


0.95% 


0.46% 


0.22% 
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Both the spectral method and the Monte Carlo algorithm have been implemented in VB.NET 
on a Pentium M processor with 1 GB of memory. 

The spectral and Monte Carlo methods agree well on the prices of variance swaps because the 
corresponding payoff is linear in the realized variance. This is to be expected since the construction 
in Section [2] of the continuous-time chain If implies that the infinitesimal drifts of the processes 
If and Tif coincide (see ([I]) and ([3])), which means that the random variables It and 'St are equal 
in expectation. The discrepancy between the two methods in the case of non-linear payoffs arises 
because the higher infinitesimal moments of It and do not coincide. 

Note also that, as mentioned in the introduction, the spectral method allows the calculation 
of the Greeks (delta, gamma, vega) of the volatility derivatives without adding computational 
complexity. In other words the gamma of the variance swap, for example, can be obtained for any 
of the maturities above in approximately 130 seconds using the same simple algorithm as in 15.11 
On the other hand it is not immediately clear how to extend the Monte Carlo method to find the 
Greeks without substantially adding computational complexity. 

6. Conclusion 

In this paper we introduce an approach for pricing derivatives that depend on pure realized 
variance (such as volatility swaps and variance swaptions) and derivatives that are sensitive to the 
implied volatility smile (such as forward-starting options) within the same framework. 

The underlying model is a stochastic volatility model with jumps that has the ability to switch 
between different CEV regimes and therefore exhibit different characteristics in different market 
scenarios. The structure of the model allows it to be calibrated to the implied volatility surface 
with minimal explicit time-dependence. The stationary nature of the model is best described by 
the implied forward smile behaviour that can be observed in Figures El El [9] and [TOl 

The model is then extended in such a way that it captures the realized variance of the underlying 
process, while retaining complete numerical solubility. Two key ideas that make it possible to keep 
track of the path information numerically are: 

• the observation that path-dependence can be expressed as the lifting of a Markov generator, 
and that 

• the state-space of the lifted process can be chosen to be on the circle so that numerical 
tractability is retained. 

Having obtained the joint probability distribution function for the realized variance and the un- 
derlying process, we outline the pricing algorithms for derivatives that are sensitive to the realized 
variance and the implied volatility within the same model. 

Appendix A. Diagonalization algorithm for partial-circulant matrices 

In this section our aim is to generalize a known diagonalization method from linear algebra which 
will yield a numerically efficient algorithm for obtaining the joint probability distribution of the 
spot and realized variance at any maturity. Let us start with some well-known concepts. 
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A matrix C G 
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C2 
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Cl 
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where each row is a cychc permutation of the row above. The structure of matrix C can also be 
expressed as 

where Cij is the entry in the i-th row and j-th cohimn of matrix C. It is clear that any circulant 
matrix is a Toeplitqj operator and in fact circulant matrices are used to approximate general 
Toeplitz matrices and explain the asymptotic behaviour of their spectra. We will not investigate 
this idea any further (for more information on the topic see [3]) since our main interest lies in a 
different generalization of circulant matrices, namely that of partial-circulant matrices, which will 
be defined in subsection lA.2[ Before doing that we are going to recall some of the known properties 
of circulant matrices. 



A.l. Eigenvalues and eigenvectors of circulant matrices. Let C denote a circulant matrix 
of dimension n as defined above. The eigenvalue A and the eigenvector y £ M" are the solutions 
of the equation Cy = Xy, which is equivalent to the system of n linear difference equations with 
constant coefficients: 

n-l 

'^Ckyk = At/o and 

k=o 

i-1 n-l 

^Cn-j+kVk + '^Ck-jyk = Xyj, for j G {1, . . . , n - 1}. 

k=Q k=j 

The variables yk, for k G {0, ...,n — 1}, in these equations are simply the coordinates of the 
eigenvector y. Such systems are routinely solved by guessing the solution and proving that it is 
correct (see appendix 1 in [20]). The solution in this case is of the form 



A = CkZ and yj = 

k=0 




for j G {0,...,n-l}, 



where z is a complex number which satisfies = 1. This implies that the eigenvalue-eigenvector 
pairs of matrix C are parameterized by the n-th roots of unity which are of the form — 
exp(— 27rir/n), where the index r lies in {0,... ,n — 1} and i is the imaginary unit. Therefore 
the j-th. coordinate of the r-th eigenvector, together with the corresponding eigenvalue, can be 



For definition see for example [1] . Toeplitz operators arise in many contexts in tlieory and practice and tlierefore 
constitute one of tfie most important classes of non-self-adjoint operators. They provide a setting for a fruitful 
interplay between operator theory, complex analysis and Banach algebras. 
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expressed as 



(7) y'f^ = -^e~'^'^ and 



1 



3 



/n 

n-l 



(8) A, = ^^Cfce-^^^'^ for r, j G {0, . . . , n - 1}. 

fc=0 

This representation is extremely useful because it allows us to deduce a number of fundamental 
facts about circulant matrices. Let us start with the eigenvectors. It is obvious that if we put all 
n vectors y^'''^ side by side into a matrix, the determinant of the linear operator obtained is the 
Vandermonde determinant, which is non-zero since its parameters are the n distinct solutions of 
the equation = 1. This implies that matrix C can be diagonalized and that all its eigenvectors 
are of the form ([7]). 

Another key property of circulant matrices is that they can all be diagonalized using the same 
set of eigenvectors. This follows directly from ([7]) since the expression for the vectors y^"^^ are clearly 
independent of matrix C . 

Expression ([8]) tells us that the r-th eigenvalue of C equals the value (at the point r) of the 
discrete Fourier transform (DFT) of the sequence (cj)j=o,...,n-i- We can therefore recover the 
sequence (cj) from the spectrum {\r)r=o,...,n-i of C using the inverse discrete Fourier transform. 
Even though this is a very well-known and celebrated fact, we will now present a short proof for 
it, because the argument itself sheds light on the behaviour of circulant matrices. 

Note that for any index G Z, such that [k mod n) is different from zero, we obtain 

r=0 1 - e n 

by summing a finite geometric series. In particular this implies that the above sum for any /c G Z 
equals ndi^k mod where 6 is the Kronecker delta which takes value 1 at zero and value everywhere 
else. The inversion formula for the DFT is now an easy consequence 

^ n— 1 ^ n— In— 1 

r=0 r=0 k=0 

A;=0 r=0 

for any I in {0, . . . ,n — 1}. Before proceeding we should note that the argument we have just 
outlined implies that a circulant matrix is uniquelj0 determined by its spectrum. 

Another consequence of the extraordinary identity Q is that for any pair of distinct indices k 
and r in {0, . . . , n — 1}, the corresponding eigenvectors y^''^ and y^'-* are perpendicular to each other. 
Since we have chosen the vectors y^^^ in ([7|) so that their norm is one, the set of all eigenvectors of 
a circulant matrix is an orthonormal basis of the vector space C"'. 

Let A be another circulant matrix given by the sequence {ak)k=o,...,n-i with the spectrum 
{ar)r=o,...,n-i- Since A and C can be diagonalized simultaneously using the basis {y^^'^; r = 
0, . . . , n — 1}, it follows that the product AC is also diagonal in this basis and that its eigenvalues 
are of the form Oj-Xr- Therefore AC is a circulant matrix whose first row equals the convolutior 



Such a statement is untrue even for self-adjoint and unitary operators. 
^Recall that the DFT of the convolution of two sequences equals the product of the DFTs of each of the sequences. 
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of the sequences (a^) and (c^). The diagonal representation also implies that the matrices A and 
C commute. Finally note that the sum ^ + C is also a circulant matrix. 

A. 2. Partial-circulant matrices. We are now going to define a class of matrices, that will include 
the Markov generator given by which can be diagonalized by the semi-analytic algorithm from 
subsection IA.4I 

Let A be a linear operator represented by a matrix in and let B^^\ for A; = 0, . . . , m — 1, 

be a family of n-dimensional matrices with the following property: there exists an invertible matrix 
U G C'^^'^ such that 

f/-i5(fe)f/ ^ ^(fc)^ for all A: G {0,...,m- 1}, 

where A^*^^ is a diagonal matrix in C"^". In other words this condition stipulates that the family of 
matrices B^^^ can be simultaneously diagonalized by the transformation U. Therefore the columns 
of matrix U are eigenvectors of B^^"^ for all k between and m — 1. 

Let us now define a large linear operator A, acting on a vector space of dimension mn, in the 
following way. Clearly matrix A can be decomposed naturally into blocks of size n x n. Let 
Aij denote an n x n matrix which represents the block in the i-th row and j-th. column of this 
decomposition. We now define the operator A as 

(10) Au := B^'^+AulRn and 

(11) Aij := Aijl^n, for all i, j £ {1, . . . ,m} such that i ^ j- 

The real numbers Aij are the entries of matrix A and Ik™ is the identity operator on M". We may 
now state our main definition. 

Definition. A matrix is termed partial-circulant if it admits a structural decomposition as in (jlOp 
and (jlip for any matrix A G R™^*" and a family of n-dimensional circulant matrices B^''\ for 
k = 0, . . . ,m — 1. 

The concept of a partial-circulant matrix is well-defined because, as we have seen in subsec- 
tion IA.lt any family of circulant matrices can be diagonalized by a unitary transformation whose 
columns consist of vectors y^^\ for r between and n — 1 (see equation ([7])). 

The operator A is very big indeed. The typical values that are of interest to us for the dimensions 
m and n are 150 and 200 respectively. This implies that matrix A contains (150-200)^, i.e. around 
one billion, entries. This means that even storing A on a computer requires about 10 Gb of memory. 

Our task is to find the spectrum of the operator A. Given its size and the fact that it is 
not a sparse matrix, this problem at first sight appears not to be tractable. But the structure 
of matrix A, combined with the ubiquitous idea of invariant subspaces of linear operators, will 
yield the solution. We will describe the diagonalization algorithm for partial-circulant matrices in 
subsection IA.4I Before we do this we need to recall the basic properties of invariant subspaces. 

A. 3. Invariant subspaces of linear operators. Let A : V ^ V he a linear operator on a 
finite-dimensional vector space V. By definition a subspace X of y is an invariant subspace of the 
operator A if and only if AX C X. Note that the set AX is a subspace V. It is clear from the 
definition that vector spaces {0}, V, AV and ker(^) are all invariant subspaces of the operator A. 
Another trivial example is the space of all eigenvectors of A that belong to an eigenvalue A. 
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It is the non-trivial examples however that make this concept so powerful. If we can find two 
invariant subspaces Xi and X2 of V for the operator A, such that Xi D X2 = {0} and dim Xi + 
dimX2 = dimy (i.e. V = Xi © X2), then in the appropriate basis the matrix representing the 
operator A takes the form 

where Ai (resp. A2) is the matrix acting on the subspace Xi (resp. X2). The zeros in the 
above expression represent trivial linear operators that map the subspace Xi into the origin of the 
subspace X2 and vice versa. 

The advantage of this structural decomposition of the original operator A is clear because it 
reduces the dimensionality of the problem. The spectral decomposition (i.e. the eigenvalues and 
eigenvectors) of A can now be obtained from the spectral decomposition of two smaller operators 
Ai and A2- Block- diagonalization consists of finding the transition matrix F (i.e. the appropriate 
coordinate change) that will transform the original matrix A into block-diagonal form given by 
matrix D above: 

P-^AF = D. 

A. 4. Algorithm for block-diagonalization. Let A be the linear operator defined in (fTUl) and ([TT]) 

which acts on the vector space C"^". We are now going to describe the block-diagonalization algo- 
rithm for A. In other words we are going to find invariant subspaces Vj of the operator A (where 
j ranges between 1 and n), such that C"*^ = Vi © • • • © Vn, and a transition matrix F e C""^^'""', 
such that the only non-zero n x n blocks of matrix F~^AF are the diagonal ones. 

Recall that, by definition of A, there exists a matrix U G (C"-x" consisting of eigenvectors for the 
matrices B^^\ Put differently the columns uj G C", for j = 1, . . . , n, of U satisfy the identity 

= Xf\j for k£{0,...,m- 1}. 

(7) 

Now fix any index j between 1 and n and define vectors £ C™"", where i = 1, . . . , m, as follows: 
(12) vl'^ := {0^_^,u'j,0^_^y, 

(i— l)n (m— i)n 



U 



'J- 



where u'j is a row of n complex numbers obtained by transposing and conjugating the vector 

(7) 

We can now define the subspace Vj of C™"" as the linear span of vectors . It is clear that the 
intersection of subspaces Vj and is trivial for any two distinct indices j,k £ {1, . . . ,n}. This is 
because the eigenvectors uj and are linearly independent in C" since, by assumption, matrix U 
is invertible. It follows directly from the definition that the dimension of Vj is m. Since there are 
exactly n subspaces Vj, we obtain the decomposition C™" = Fi © • • • © Ki- 

If we manage to show that each space Vj is an invariant subspace for the operator A, we will be 
able to conclude that A can be expressed in the block-diagonal form as described in subsection lA. 31 
Since the subspace Vj is defined as a linear span of a set of vectors {vp^ i = l,...,m}, the 
invariance property AVj C Vj will follow if we demonstrate that the vector Av^''^ is in Vj for all 
i = 1, . . . ,m. By definition of A ( (|10p and (jlip ) it immediately follows that 

m 

(13) Avl'^ =Y,A^kvl;!^ + for ie{l,...,m}, 

k=l 
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where A^* is the eigenvalue of matrix B^^~^^ that corresponds to the eigenvector n,- and the real 
numbers A^i are the entries of matrix A G R™^™. Identity (I13p implies that each subspace Vj is 
an invariant subspace for A. Furthermore, if we define a matrix F G ([^"'■nxmn -j^ ^j^g following way 

(14) F:=(.f\...,.«,.f\...,.W,...,.i"\...,.W), 

then matrix D = F~^AF is block-diagonal. In other words if we decompose D into matrices 
Dij of size m X m, then the following formula holds 

(15) Aj =5.,(A + e(j)) for i,jG {!,..., n}, 

where Q^^^ is a diagonal matrix in C'"^™ with its i-th diagonal element equal to A^* ^\ As usual 
the symbol 5ij denotes the Kronecker delta function. 

Expression p5|) gives us the block-diagonal representation of the operator A. Notice that the 
diagonal elements of matrix Q^^^ are precisely the eigenvalues of matrices i?*-*-*, for i = 0, . . . , m — 1, 
that correspond to the eigenvector Uj. 

The algorithm to block-diagonalize the operator A, defined by matrices A G R"^^™ and B^^'^ G 
j^nxn |-ggg pQj) and ([TT]) ). can now be described as follows: 

(I) Find matrix U G C"^" whose columns are the common eigenvectors uj, for j G {1, . . . ,n}, 
of the family B^''\ 

(II) Construct the transition matrix F using the columns of matrix U as described in (1120 
and (fTil) . 

(III) Find the eigenvalues A^'^^ which satisfy B^^^uj = ^j^^Uj for k G {0, ...,m — 1} and j G 
{l,...,n}. 

(IV) Construct diagonal matrices Q'^^^ G C"*^™, for all j G {1, . . . ,n}, given by B^-^^ = ^ifcA^* ^\ 
where the indices i, A; run over the set {!,..., m}. 

(V) Construct the block-diagonal representative D for the operator A as described in (|15p . 

Our main task is to find the spectrum of matrix A. Notice that, since the spectrum of A is 
a union of the spectra oi A-\- this algorithm has reduced the problem of diagonalizing an 
nm X nm matrix to finding the spectra of n matrices of size m x m. The algorithm provides a key 
step for our pricing method because it enables us to model the behaviour of the realized variance 
by increasing the numerical complexity only linearly . 

We should also note that in the case of the lifted Markov generator L in matrix A is the 
generator L of the underlying process while the family B^'^^ consists of circulant matrices. In other 
words the operator L is given by a partial-circulant matrix. It therefore follows from the discussion 
in subsection I A . 1 1 that the columns of the corresponding transition matrix F are pairwise orthogonal 
and that the entries of matrix LF are the values of a partial discrete Fourier transform of the rows 
of L. This simple observation is useful when calculating the probability kernel of the lifted generator 
in Section [3l 



Since each row of matrix L is naturally described by two variables, namely the value of the underlying and the 
value of the realized variance, partial DFT is by definition a DFT acting on the second variable. 
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A. 5. Proof of Theorem 13.11 Let us start by recalling that any holomorphic function defined on 
C has a Taylor expansion around zero that converges everywhere. We can therefore define (/'(^) 
using the power series of the function (p for any linear operator A on a finite-dimensional vector 
space. It also follows from the fact that (j) has a Taylor expansion that any invariant subspace (see 
subsection IA.3I for definition) of A is also an invariant subspace of (p{A). In particular if A has 
a block-diagonal decomposition in the sense of appendix \AA\ then the matrix (p{A) also has one. 
Moreover if i3 is a block in A, then (l){B) must be a block in (j){A). 

We know that the Markov generator L can be expressed as L = FDF^^, where Z) is a block- 
diagonal matrix of the form 

/Lo ••• \ 
Li ••• 

D = 

^0 ••• L2cJ 

and the transition matrix F is given by (I14p . We have just seen that (j){D) must therefore also be 
in block-diagonal form: 



( '/'(Lo) 






'/'(Li) 








\ 



Y ••• (PiUc)) 
The power series expansion of 4> implies that 0(L) = F(j){D)F~^ . Since matrix F is defined using 
the eigenvectors of circulant matrices, it follows immediately that the inverse F~^ can be obtained 
by transposing F and conjugating each of its elements. Note that the dimension of our circulant 
matrices is 2C + 1 and express (p(L){x, P,c;y,^,d) = {u,(j){D)v) as a real inner product of two 
vectors u and (p{D)v, where u equals the (rr, /3, c)-row of the matrix F and v is the (y, 7, (i)-column 
of (i.e. the conjugated (y,7,(i)-row of F). 

It follows from the definition of F and the above expression for (f>{D) that the non-zero coordinates 
of the vector u are of the form 

V2C + 1 

for all /c G ^ and that the corresponding coordinates of (j){D)v are 

The equality in the theorem now follows directly from the expressions for the coordinates of the 
vectors u and <j){D)v and the definition of the real inner product. This concludes the proof of 
Theorem 13.11 

Appendix B. Volatility derivatives 

In this appendix we are going to give a brief description of the volatility derivatives discussed 
in this paper. We start with the simplest case, namely a forward on the realized variance, which 
defines a variance swap. In subsection IB . 2 1 we define options with payoffs that are general functions 
of realized variance. Subsection IB. 31 concerns derivatives that are dependent on implied volatility. 
In particular we recall the definition of the forward-starting options and of the implied volatility 
index. 



SPECTRAL METHODS FOR VOLATILITY DERIVATIVES 



25 



B.l. Variance swaps. As mentioned above, a variance swap expiring at time T is simply a forward 
contract on the realized variance Sy, quoted in annual terms, of the underlying stock (or index) 
over the time interval [0,T]. The payoff is therefore of the form 

(St - Kvar)A^, 

where -ftTvar is the strike and N is the notional of the contract. The fair value of the variance is the 
delivery price i^Tvar which makes the swap have zero value at inception. 

At present such contracts are liquidly traded for most major indices. The delivery price is usually 
quoted in the markets as the square of the realized volatility, i.e. -f^var = where i^T is a value of 
realized volatility expressed in percent. The notional N is usually quoted in dollars per square of 
the volatility poin 

10. 

A key part of the specification of a variance swap contract is how one measures the realized 
variance S^- There are a number of ways in which discretely sampled returns of an index (or of 
an index futurj^ Ft) can be calculated and used for defining the realized variance. We will now 
describe the two most common approaches. 

The usual definition of the annualized realized (i.e. accrued) variance of the underlying process Ft 
in the period [0, T], using logarithmic returns, is - X^"=i(log J^*^ )^, where times tj, for i = 0, . . . , n, 
are business days from now to = until expiry t„, = T. The normalization constant d is the number 
of trading days per year. Another frequently used definition of the realized variance is given by 
nZ]r=i(~T — ^^)^- It is a standard fact about continuous square-integrable martingales that in 
the limit, as we make partitions of the interval [0, T] finer and finer, both sums exhibit the following 
behaviouiF^: 

(logF), = hm flog^)' = hm f ^V^'^O ^- 

The convergence here is in probability and the process (log-F)j is the quadratic- variatiorF^ process 
associated to log(Ft). 

In our framework the underlying process Ft is a continuous-time Markov chain (see chapter 6 
of [27j for definitions and basic properties). We define the annualized realized variance (of Ft 
over the time interval [0,T]) to be the limit 

(16) ^t:=1 hm yf^^^V^'' 

^ ' TpH^o^^K Ft,_, 



10a 

volatility point is one basis point of volatility, i.e. 0.01 if volatility is quoted in percent. This means that the 
quote for the notional value of the variance swap tells us how much the swap owner gains if the realized variance 
increases by 0.0001 = 0.01^ 

^^The reason for considering index futures rather than the index itself is twofold. The futures are used for hedging 
options on the index because they are much easier to trade than the whole portfolio of stocks that the index comprises. 
Also, it is well-known that futures prices are martingales under the appropriate risk-neutral measure which depends 
on the frequency of mark-to-market. If the futures contract marks to market continuously, then the price process Ft 
is a martingale in the risk-neutral measure induced by the money market account as a numeraire. Otherwise we have 
to take the rollover strategy, with the same frequency as mark-to-niarket, as our numeraire to obtain the martingale 
measure for Ft. 

Notice also that these equalities hold because the difference of the process log(Ft) and Jl ^ is of finite 
variation, which is a consequence of Ito's lemma (see Theorem 3.3 in [24]). 

^■^For a precise definition of a quadratic- variation of continuous square-integrable martingales see chapter 1 of [24| . 
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where, for every n G N, the set (tj)i=o,...,n is a strictly increasing sequence of times between and 
T and p{n) := maxjtj — i = 1, . . . ,n} is the size of the maximal subinterval given by the 
sequence (ti)i=o,...,n (cf- Definitions ([T]) in Section [5] and ^ in Subsection 12. It should be 
noted that the techniques described in Sections [2] and [Aj which provide numerical solubility for 
our model, can be generalized to the situation where the realized variance is defined as a discretely 
sampled sum in ()16p but without the limit. We are not going to pursue this line of thought any 
further, but should notice that the discrete definition of the realized variance would require the 
application of the block-diagonalization algorithm (appendix [X]) to the probability kernel between 
any two consecutive observation times ti rather than the application of the algorithm to the Markov 
generator directly, which is what is done in Section O 

B.2. General payoffs of the realized variance. A volatility swap is a derivative given by the 
payoff 



where aj! is the realized volatility over the time interval [0, T] quoted in annual terms, -ftTvoi is the 
annualized volatility strike and N is the notional in dollars per volatility point. The market con- 
vention for calculating the annualized realized volatility cj|? differs slightly from the usual statistical 
measur of a standard deviation of any discrete sample and is given by the formula 



where d is the number of trading days per year and ti are business days from now to = until 
expiry t„ = T of the contract. For our purposes we shall define realized volatility Uj;, quoted in 
annual terms, over the time interval [0, T] as 



where St is the annualized realized variance defined in (|16|) . It is clear from this definition that 
the payoff of the volatility swap can be view as a non-linear function of the realized variance. 

Since volatility swaps are always entered into at equilibrium, an important issue is the determi- 
nation of the fair strike K^oi for any given maturity T. As discussed in Section [H a term structure 
of such strikes must be part of the market data that some models require (e.g. [7]) in order to be 
calibrated. In our case the strikes K^^i, for any maturity, are implied by the model which uses as 
its calibration data the market implied vanilla surface. The value of K^^i for a given maturity T 
is then given by the expectation Eo[\/Sr], which can easily be obtained as soon as we have the 
probability distribution function for T,t- 

The same reasoning applies to variance swaps. The fair strike i^var for the variance swap of 
maturity T can, within our framework, be obtained by taking the expectation Eo[Sj']. It therefore 

^^Given a sample of n values Xi, . . . , X„ with the mean A* = ^ "127=1 ^^e unbiased statistical estimation of the 
standard deviation is given by 




C^T '■= 



n 
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follows from the concavity of the square root function and Jensen's inequalit}0 that the following 
relationship holds between the fair strikes of the variance and volatility swaps 

for any maturity T. This inequality is always satisfied by the market quoted prices for variance 
and volatility swaps and is there to account for the fact that variance is a convex function of 
volatility. Put differently, this is just a convexity effect, similar to the one observed for ordinary 
call options, related to the magnitude of volatility of volatility. The larger the "vol of vol" is, the 
greater the convexity effect becomes. This phenomenon can be observed clearly in the markets with 
a very steep skew for implied volatilities. If one wanted to estimate its size in general, it would be 
necessary to make assumptions about both the level and volatility of the future realized volatility. 
Within our model this can be achieved directly by comparing the values of the two expectations 
(see Figure (llSp for this comparison based on the market implied vanilla surface for the S&P 500). 

There are other variance payoffs which are of practical interest and can be priced and hedged 
within our framework. Examples are volatility and variance swaptions whose payoffs are (-^/S^^ — 
-^voi)^ and — Kvar)"*" respectively, where as usual (2;)+ equals max(2;,0) for any 2; G M. Capped 
volatility swaps are also traded in the markets. Their payoff function is of the form (min(-v/Sr, am) — 
Kvoi), where am denotes the maximum allowed realized variance. It is clear that all such contracts 
can be priced easily within our framework by integrating any of these payoffs against the probability 
distribution function (see Figure [TT] for maturities above 6 months) of the annualized realized 
variance T,t and then multiplying the expectation with the corresponding discount factor. 

It should be noted that there exist even more exotic products, like corridor variance swaps 
(see [TU] ) , whose payoffs depend on the variance that accrues only if the underlying is in a predefined 
range. Such products cannot be priced directly in the existing framework. A modification of the 
model would be required to deal with this class of derivatives. However we will not pursue this 
avenue any further. 

B.3. Forward staring options and the volatility index. Let T' and T be a pair of maturities 
such that T' < T. A forward staring option (or a forward- start) is a vanilla option with expiry T 
and the strike, set at time T' , which is equal to aSx'- The quantity St is the underlying financial 
instrument the option is written on (usually a stock or an index). More formally the value of a 
forward-start at time T (i.e. its payoff) is given by 

(17) Vfs{T) = {St - aST')^ , 

where the constant a is specified at the inception of the contract and is know as the forward strike. 
It is clear from the definition of the forward-start that its value at time T' equals the value of a 
plain vanilla call option 

VFs{T') = Vc{ST',T-T',aST') 
that expires in T — T' years and whose strike equals aSx'- Notice that at time T' the constant 
a can be characterized as the ratic0 between the spot price St' and the strike of the call option 
into which the forward-start is transformed. In the classical Black-Scholes framework, we have an 

"'^^For any convex function </) : R — > R and any random variable X : SI R with a finite first moment Jensen's 
inequality states that <;/)(E[X]) <,¥.[(j){X)\. 

-'^''This ratio is sometimes referred to as the "moneyness" of the option. 
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explicit formula, denoted by BS{St',T — T' , aST',r', a'), for the value of this call option (see [3]). 
This formula depends linearly on the spot level St' if the ratio of the spot and the strike (i.e. 
the "moneyness") is known. In other words, assuming we are in the Black-Scholes world with a 
deterministic term-structure of volatility and interest rates and zero dividends, we can express the 
value of the forward-start at time T' as 

Vfs{T') = St'BS{1, T - T', a, r', a'), 

where a' is the forward volatility rat and r' is the forward interest rate between T' and T. The 
following key observations about the Black-Scholes value and sensitivities of the forward-starting 
option are now clear: 

• the value equals Vfs(O) = 5oBS(l, T — T' , a, r' , a') and the delta (i.e. -^^Fsi^)) is simply 
BS{1,T -T',a,r',a'), 

• the forward-starting option is gamma neutra (i.e. ^Vfs{0) = 0) and 

• the contract has non-zero vega (i.e. ^Vfs(O) > 0). 

We are interested in the Black-Scholes pricing formula for the forward-starts because we need to 
use it when expressing the forward volatility smile of our model. The values of forward volatility 
a', implied by the equation Vfs{0) = SoBS{l,T — T' , a,r' , a'), are plotted in Figures [71 [HI [9l and [TOl 
for a wide range of values of the forward strike a and a variety of time horizons T', T. In other 
words, we first calculate the value Vpsi^) of the forward-start and then invert the Black-Scholes 
pricing formula to obtain the implied forward volatility a' . 

We shall now give a brief description of the implied volatility index (VIX) and then move on 
to discuss the future probability distribution of VIX, which, as will be seen, is directly related to 
forward-starting options. 

VIX was originally introduced in 1993 by Chicago Board Options Exchange (CBOE) as an index 
reflecting the 1 month implied volatility of the at-the-money put and call options on S&P 100. To 
facilitate trading in VIX, in 2003 CBOE introduced a new calculation, using a full range of strikes 
for out-of-the-money options, to define the value of VIX. At the same time the underlying financial 
instrument on which the options are written was changed to S&P 500 (for a detailed description of 
these changes and their ramifications see |12j). The new formula is 

(18) 4„ = |^Mi,.T<j,;,,)_|(^£_iy_ 

where the index itself is given by VIX = lOOcvix- The sequence Ki consists of all the exchange 
quoted strikes and the quantities Q{Ki) are the corresponding values of out-of-the-money put/call 
options expiring at maturity T, where T equals 1 monthThe quantity F in formula ()18p is the 
forward value of the S&P 500 index derived from option price^lj and the at-the-money strike Kq 
is defined as the largest strike below F. Note also that it is precisely at Kq that the symbol Q{K) 
in (jlSp changes from put to call options. 

-'^ ''^Assuming that the term-structure of volatility is parametrized by (j{t), the forward volatility rate a' is given by 

^^It should be noted that, in the presence of stochastic volatility, the gamma of a forward-starting option is no 
longer necessarily zero. However in a realistic model it should not be too large because it reflects the dependence of 
volatility on very small moves of the underlying, the effect of which should be negligible. 

^^See page 3 in [T^ for the precise definition. 
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The reason why formula (jlSp ahows easier trading of volatihty follows from the simple obser- 
vation that CTyj-j^ is essentially the value of a European derivative, expiring at time T, with the 
logarithmic payoff given in (|19|) . This is a consequence of the well-known decomposition of any 
twice differentiable payoff described in [5], [H], [T^ and other sources: 

(19) - log = -^^y^ + lJ^iK- Sr^dK + -^(5^ - K)+dK. 

This formula holds for any value of but expressions simplify if we assume that F equals the 
forward of the index St at time T (i.e. F = 'K[St\)- By taking the expectation with respect to the 
risk-neutral measure we get the following expression for the forward price of the log payoff: 

(20) - E [log = ±e^Tp{K)dK + ^e'-^C7(if)dif, 

where C{K) = e-^^E[(5T - K)+] (resp. P{K) = e-^^E[(K - 5t)+]) is the price of a call (resp. 
put) option struck at K. It is shown in [15] that the portfolio of vanilla options given by (j20p 
can be used to hedge perfectly a variance swap if there are no jumps in the underlying market. 
From our point of view the expression (I20p is interesting because a simple calculation shows that 
definition ()18|) is a possible discretization of it. By defining its own version of the approximation to 
the logarithm, CBOE has created a volatility index which can be replicated by trading a relatively 
simple European payoff. This feature greatly simplifies the trading of VIX. 

Since the implied volatility index is defined by a portfolio of puts and calls in (fTSj) . it is clear 
that the random nature of the value of VIX at time t will be determined by the value of the 
corresponding portfolio of forward-starting options. The probability distribution function for the 
behaviour of the volatility index at t is obtained from the model by the following procedure: 
(I) Fix a level S of the underlying. 

(II) Find the probability that the price process is at level 5" at a given time t in the future. 

(III) Evaluate the portfolio of options that define the volatility index between times t and t + T, 
conditional on the process St being at level S. 

(IV) Repeat these steps for all attainable levels S for the underlying Markov chain St- 

(V) Subdivide the real line into intervals with disjoint interiors of length 5, where 5 is a small 
positive number. To each of the intervals assign a probability that is a sum of probabilities 
in step (II) corresponding to the values obtained in step (III) that lie within the interval. 

This describes the construction of the probability distribution function of the volatility index at 
time t. The plot of this pdf for a variety of maturities can be seen in Figure [T9l 

Our final task is to price any European payoff written on the level of VIX at a certain time 
horizon. Given that our model allows us to extract the pdf of VIX for any expiry, pricing such a 
derivative amounts to integrating the payoff function against the probability distribution that was 
described above. 
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Figure 1. implied volatilities for out-of-the-money European options written on the S&P 500 
with maturities between 3 months and 10 years. For each market-specified maturity the squares 
■ represent market-implied volatilities. The continuous curves graph the implied volatility of the 
model as a function of strike. The relative value of S&P 500, with respect to the current level of 
spot, is plotted along the line of abscisse. 
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Figure 2. implied probability density function for the S&P 500 under the forward measure. 
The value of the S&P 500 plotted along the line of abscisse is a relative value (in percent) with 
respect to the current level of spot. The pdfs between 6 months and 5 years can also be viewed as 
the rescaled marginals of the joint pdfs in figures 1151 1161 and 1171 
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Calendar time 

Figure 3. Deterministic time-change f(t) (measured in years) as a function of calendar time t 
(also in years). Function / was used in the cahbration of the model to the vanilla surface of the 

S&P 500. 
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Figure 4. Delta profiles of call options on the S&P 500 with maturities between 6 months and 
2 years, all struck at 100. We calculate the delta of a call option (i.e. A(S) = §§(5')) for all lattice 
points S using a symmetric difi'erence as described in Subsection [STT] Notice that we are using the 
underlying index 5 and the strike of the options on their relative scales with respect to the level at 
which the index was trading when the snapshot of the market was taken. 
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Figure 5. Gamma profiles of call options on the S&P 500 with maturities between 6 months 
and 2 years, all struck at 100. We calculate the gamma of a call option (i.e. r{S) = |^(S')) for all 
lattice points S using a symmetric difference as described in Subsection 15.11 The same comment 
as in Figure |4l about the relative value of the index and the strike, applies. 
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Figure 6. Vega profiles of call options on the S&P 500 with maturities between 6 months and 
2 years, all struck at 100. We are calculating the vega of a call option by bumping the current 
volatility regime, repricing the option and plotting the difference from the original option value, 
for all points on the lattice. Notice that vega and gamma profiles are very similar in shape but 
different in magnitude, which is consistent with the general market view on the two Greeks. 
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Figure 7. implied forward volatility skews 3 months from now. For example the code 3mly 
means that T' equals 3 months and that the expiry T can be calculated as T = (T' + 1 year) , where 
T' ,T are as defined in Subsection lB.3l Along the line of abscisse we plot the forward strike a (i.e. 
the "moneyness" of the ordinary call option that the forward-start becomes at time T') as defined 
in formula (fT7)l . The ordinate axis contains the forward volatility values expressed in percentage, 
as impled by the model. 
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Figure 8. implied forward volatility skews 6 months from now. For example the code 6mly 
means that T' equals 6 months and that the expiry T can be calculated as T = (T' + 1 year) , where 
r',r are as defined in Subsection lB.3l Along the line of abscisse we plot the forward strike a (i.e. 
the "moneyness" of the ordinary call option that the forward-start becomes at time T') as defined 
in formula ([17)1 . The ordinate axis contains the forward volatility values expressed in percentage, 
as impled by the model. 
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Figure 9. implied forward volatility skews 1 year from now. For example the code ly6m means 
that T' equals 1 year and that the expiry T can be calculated as T = (T' + 6 months), where T' ,T 
are as defined in Subsection IB. 31 Along the line of abscisse we plot the forward strike a (i.e. the 
"moneyness" of the ordinary call option that the forward-start becomes at time T') as defined in 
formula (|17|l . The ordinate axis contains the forward volatility values expressed in percentage, as 
impled by the model. 
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Figure 10. implied forward volatility skews 2 years from now. For example the code 2y6m 
means that T' equals 2 years and that the expiry T can be calculated as T = (T' + 6 months) , 
where T',T are as defined in Subsection IB. 31 Along the line of abscisse we plot the forward strike 
a (i.e. the "moneyness" of the ordinary call option that the forward-start becomes at time T') 
as defined in formula p7p . The ordinate axis contains the forward volatility values expressed in 
percentage, as impled by the model. 
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Figure 11. Probability distribution functions for the realized variance, quoted in annual terms, 
for maturities between 6 months and 5 years as given by our model after calibration to the implied 
volatility surface of the S&P 500 (see Section |4] for details). These pdfs are marginal distributions 
obtained from the joint probability distribution function ([6]) by integrating it in the dimension of 
the spot value of the index. 



Probability distribution for the realized variance in a stochastic volatility model without jumps 
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Figure 12. Probability distribution functions for the realized variance, quoted in annual terms, 
for maturities between 6 months and 5 years as given by the model with zero jump intensities, 
based on two volatility regimes. For the full list of values of the parameters used to specify this 
model see Subsection 15.51 
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Figure 13. Probability distribution functions for the realized variance, quoted in annual terms, 
for maturities between 6 months and 5 years as given by the model used in Figure 1111 but with 
zero jump intensities, based on two volatility regimes. The model parameters are the same as the 
ones given in table [T] of Section but with — 0, for a = 0, 1. 
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Figure 14. The spacing of the grid a{t) = f{t) for the realized variance as a function of 
calendar time t. The constant C equals 100, for all times t from now to 5 years, and corresponds 
to 201 states for the process It that represents the realized variance of the underlying up to time t. 
The function f{t) represents financial time and is graphed in Figure |3] 
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Figure 15. Joint probability distribution function for the annualized realized variance and the 
spot rate of S&P 500 in 6 months' and 1 year's time. 
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Figure 16. Joint probability distribution function for the annualized realized volatility and the 
spot rate of S&P 500 in 2 years' and 3 years' time. 




Figure 17. Joint probability distribution function for the annualized realized volatility and the 
spot rate of S&P 500 in 4 years' and 5 years' time. 
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Term structures of variance swap prices and related instruments 
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Figure 18. Term structures of variance swap prices for maturities between 6 months and 5 years 
as implied by the vanilla market data using the calibrated model (for precise values see table dj. 
The delivery prices (i.e. fair strikes) are computed as Eo[St], where Et is the annualized realized 
variance for each tenor T. Everything is expressed in terms of volatility, i.e. the prices are in 
percent and are obtained by taking the square root of the variance. We also plot fair delivery prices 
for volatility swaps. The convexity bias implied by the model can be clearly observed. The value of 
the log contract, defined in equation (|20p . given by the underlying model is also plotted. Observe 
that this portfolio of options is always worth more than the corresponding variance swap because, 
in our model, we only allow for down jumps (see Section |4]). This behaviour is exactly as predicted 
by the analysis in (Demeterfi et al. 19996) (equation 42) when they added a single down-jump to 
the underlying process and studied its influence on the static hedge (|20|) for the variance swap. 
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Figure 19. Probability distribution functions of the implied volatility index for maturities be- 
tween 6 months and 2 years. 
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Figure 20. Probability distribution functions of a portfolio of forward-starting options as in the 
definition of VIX (see (|18|) in Subsection IB.Sp . where time t is fixed at 1 year and time T varies 
form 1 month to 2 years (see Subsection IB . 31 for the definition of a probability distribution function 
for a portfolio of forward-starts and the role of parameters t and T). Note that the lylm pdf is, 
according to our definition, the distribution of VIX in 1 year. 
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